The present invention relates to image processing systems and methods, and more particularly to image registration systems that combine two or more images into a composite image in particular the fusion of anatomical manifold based knowledge with volume imagery via large deformation mapping which supports both kinds of information simultaneously, as well as individually, and which can be implemented on a rapid convolution FFT based computer system.
Image registration involves combining two or more images, or selected points from the images, to produce a composite image containing data from each of the registered images. During registration, a transformation is computed that maps related points among the combined set of landmarks in each of the images that are to be registered. For example, when registering two MRI images of different axial slices of a human head, a physician may label points, or a contour surrounding these points, corresponding to the cerebellum in two images. The two images are then registered by relying on a known relationship among the landmarks in the two brain images. The mathematics underlying this registration process is known as small deformation multi-target registration.
In the previous example of two brain images being registered, using a purely operator-driven approach, a set of N landmarks identified by the physician, represented by x, where i=1 . . . N, are defined within the two brain coordinate systems. A mapping relationship, mapping the N points selected in one image to the corresponding N points in the other image, is defined by the equation u(xi)=ki, where i=1 . . . N. Each of the coefficients, ki, is assumed known.
The mapping relationship u(x) is extended from the set of N landmark points to the continuum using a linear quadratic form regularization optimization of the equation:
                    u        =                  arg          ⁢                      xe2x80x83                    ⁢                                    min              u                        ⁢                          ∫                                                "LeftDoubleBracketingBar"                  Lu                  "RightDoubleBracketingBar"                                2                                                                        (        1        )            
subject to the boundary constraints u(xi)=ki,. The operator L is a linear differential operator. This linear optimization problem has a closed form solution. Selecting L=xcex1∇2+xcex2∇(∇xc2x7) gives rise to small deformation elasticity.
For a description of small deformation elasticity see S. Timoshenko, Theory of Elasticity, McGraw-Hill, 1934 (hereinafter referred to as Timoshenko) and R. L. Bisplinghoff, J. W. Marr, and T. H. H. Pian, Statistics of Deformable Solids, Dover Publications, Inc., 1965 (hereinafter referred to as Bisplinghoff). Selecting L=∇2 gives rise to a membrane or Laplacian model. Others have used this operator in their work, see e.g., Amit, U. Grenander, and M. Piccioni, xe2x80x9cStructural image restoration through deformable templates,xe2x80x9d J. American Statistical Association. 86(414):376-387, June 1991, (hereinafter referred to as Amit) and R. Szeliski, Bayesian Modeling of Uncertainty in Low-Level Vision, Kluwer Academic Publisher, Boston, 1989 (hereinafter referred to as Szeliski) (also describing a bi-harmonic approach). Selecting L=∇4 gives a spline or biharmonic registration method. For examples of applications using this operator see Grace Wahba, xe2x80x9cSpline Models for Observational Data,xe2x80x9d Regional Conference Series in Applied Mathematics. SIAM, 1990, (hereinafter referred to as Whaba) and F. L. Bookstein, The Measurement of Biological Shape and Shape Change, volume 24, Springer-Verlag: Lecture Notes in Biomathematics, New York, 1978 (hereinafter referred to as Bookstein).
The second currently-practiced technique for image registration uses the mathematics of small deformation multi-target registration and is purely image data driven. Here, volume based imagery is generated of the two targets from which a coordinate system transformation is constructed. Using this approach, a distance measure, represented by the expression D(u), represents the distance between a template T(x) and a target image S(x). The optimization equation guiding the registration of the two images using a distance measure is:
                    u        =                              arg            ⁢                          xe2x80x83                        ⁢                                          min                u                            ⁢                              ∫                                                      "LeftDoubleBracketingBar"                    Lu                    "RightDoubleBracketingBar"                                    2                                                              +                      D            ⁡                          (              u              )                                                          (        2        )            
The distance measure D(u) measuring the disparity between imagery has various forms, e.g., the Gaussian squared error distance ∫|T(h(x))xe2x88x92S(x)|2dx, a correlation distance, or a Kullback Liebler distance. Registration of the two images requires finding a mapping that minimizes this distance.
Other fusion approaches involve small deformation mapping coordinates xxcex5xcexa9 of one set of imagery to a second set of imagery. Other techniques include the mapping of predefined landmarks and imagery, both taken separately such as in the work of Bookstein, or fused as covered via the approach developed by Miller-Joshi-Christensen-Grenander the ""212 patent described in U.S. Pat. No. 6,009,212 (hereinafter xe2x80x9cthe ""212 patentxe2x80x9d) herein incorporated by reference mapping coordinates xxcex5xcexa9 of one target to a second target. The existing state of the art for small deformation matching can be stated as follows:
Small Deformation Matching: Construct h(x)=xxe2x88x92u(x) according to the minimization of the distance D(u) between the template and target imagery subject to the smoothness penalty defined by the linear differential operator L:                               h          ⁡                      (            ·            )                          =                                            ·                              -                                  u                  ⁡                                      (                    ·                    )                                                                        ⁢                          xe2x80x83                        ⁢            where            ⁢                          xe2x80x83                        ⁢                          u              ^                                =                                    arg              ⁢                                                min                  u                                ⁢                                  ∫                                                            "LeftDoubleBracketingBar"                      Lu                      "RightDoubleBracketingBar"                                        2                                                                        +                                          ∑                                  i                  =                  1                                N                            ⁢                                                D                  ⁡                                      (                    u                    )                                                  .                                                                        (        3        )            
The distance measure changes depending upon whether landmarks or imagery are being matched.
1. Landmarks alone. The first approach is purely operator driven in which a set of point landmarks xi, i=1, . . . N are defined within the two brain coordinate systems by, for example, an anatomical expert, or automated system from which the mapping is assumed known: u(xi)=ki, i=1, . . . , N. The field u(x) specifying the mapping h is extended from the set of points {xi} identified in the target to the points {yi} measured with Gaussian error co-variances xcexa3i:                               u          ^                =                              arg            ⁢                                          min                u                            ⁢                              ∫                                                      "LeftDoubleBracketingBar"                    Lu                    "RightDoubleBracketingBar"                                    2                                                              +                                    ∑                              i                =                1                            N                        ⁢                                                            (                                                            y                      i                                        -                                          x                      i                                        -                                          u                      ⁡                                              (                                                  x                          i                                                )                                                                              )                                t                            ⁢                                                ∑                                      -                    1                                                  ⁢                                                      (                                                                  y                        i                                            -                                              x                        i                                            -                                              u                        ⁡                                                  (                                                      x                            i                                                    )                                                                                      )                                    T                                                                                        (        4        )            
Here (xc2x7)r denotes transpose of the 3xc3x971 real vectors, and L is a linear differential operator giving rise to small deformation elasticity (see Timoshenko and Bisplinghoff), the membrane of Laplacian model (see Amit and Szeliski), bi-harmonic (see Szeliski), and many of the spline methods (see Wahba and Bookstein). This is a linear optimization problem with closed form solution.
2. The second approach is purely volume image data driven, in which the volume based imagery is generated of the two targets from which the coordinate system transformation is constructed. A distance measure between the two images being registered I0, I1 is defined as D(u)=∫|I0(xxe2x88x92u(x))xe2x88x92I1(x)|2dx. The corresponding optimization is:                                           h            ⁡                          (              ·              )                                =                      ·                          -                              u                ⁡                                  (                  ·                  )                                                                    ⁢                  
                ⁢        where        ⁢                  xe2x80x83                ⁢                  
                ⁢                  u          =                                    arg              ⁢                              xe2x80x83                            ⁢                                                min                  u                                ⁢                                  ∫                                                            "LeftDoubleBracketingBar"                      Lu                      "RightDoubleBracketingBar"                                        2                                                                        +                          ∫                                                                    "LeftBracketingBar"                                                                                            I                          0                                                ⁡                                                  (                                                      x                            -                                                          u                              ⁡                                                              (                                x                                )                                                                                                              )                                                                    -                                                                        I                          1                                                ⁡                                                  (                          x                          )                                                                                      "RightBracketingBar"                                    2                                ⁢                                                      ⅆ                    x                                    .                                                                                        (        5        )            
The data function D(u) measures the disparity between imagery and various forms. Other distances are used besides just the Gaussian squared error distance, including correlation distance, Kullback Liebler distance, and others.
3. The algorithm for the transformation of imagery I0 into imagery I1 has landmark and volume imagery fused in the small deformation setting as in the ""212 patent. Both sources of information are combined into the small deformation registration:                     u        =                              arg            ⁢                                          min                u                            ⁢                              ∫                                                      "LeftDoubleBracketingBar"                    Lu                    "RightDoubleBracketingBar"                                    2                                                              +                      D            ⁡                          (              u              )                                +                                    ∑                              i                =                1                            N                        ⁢                                                            "LeftDoubleBracketingBar"                                                            y                      i                                        -                                          x                      i                                        -                                          u                      ⁡                                              (                                                  x                          i                                                )                                                                              "RightDoubleBracketingBar"                                2                                            σ                1                2                                                                        (        6        )            
Although small deformation methods provide geometrically meaningful deformations under conditions where the imagery being matched are small, linear, or affine changes from one image to the other. Small deformation mapping does not allow the automatic calculation of tangents, curvature, surface areas, and geometric properties of the imagery. To illustrate the mapping problem, FIG. 9 shows an oval template image with several landmarks highlighted. FIG. 10 shows a target image that is greatly deformed from the template image. The target image is a largely deformed oval that has been twisted. FIG. 11 shows the results of image matching when the four corners are fixed, using small deformation methods based on static quadratic form regularization. These Figures illustrate the distortion which occurs with small deformation linear mapping when used with landmark points which define a motion corresponding to large deformation. As can be seen in FIG. 11, landmarks defined in the template image often map to more than one corresponding point in the target image.
Large deformation mapping produces maps for image registration in which the goal is to find the one-to-one, onto, invertible, differentiable maps h (henceforth termed diffeomorphisms) from the coordinates xxcex5xcexa9 of one target to a second target under the mapping
h:xxe2x86x92h(x)=xxe2x88x92u(x), xxcex5xcexa9xe2x80x83xe2x80x83(7) 
To accommodate very fine variations in anatomy the diffeomorphic transformations constructed are of high dimensions having, for example a dimension greater than 12 of the Affine transform up-to the order of the number of voxels in the volume. A transformation is diffeomorphic if the transformation from the template to the target is one-to-one, onto, invertible, and both the transformation and it""s inverse are differentiable. A transformation is said to be one-to-one if no two distinct points in the template are mapped to the same point in the target. A transformation is said to be onto if every point in the target is mapped from a point in the template. The importance of generating diffeomorphisms is that tangents, curvature, surface areas, and geometric properties of the imagery can be calculated automatically. FIG. 12 illustrates the image mapping illustrated in FIG. 11 using diffeomorphic transformation.
The present invention overcomes the limitations of the conventional techniques of image registration by providing a methodology which combines, or fuses, some aspects of techniques where an individual with expertise in the structure of the object represented in the images labels a set of landmarks in each image that are to be registered and techniques that use the mathematics of small deformation multi-target registration, which is purely image data driven. An embodiment consistent with the present invention uses landmark manifolds to produce a coarse registration, and subsequently incorporates image data to complete a fine registration of the template and target images.
Additional features and advantages of the invention will be set forth in the description which follows, and in part, will be apparent from the description, or may be learned by practicing the invention. The embodiments and other advantages of the invention will be realized and obtained by the method and apparatus particularly pointed out in the written description and the claims hereof as well as in the appended drawings.
To achieve these and other advantages and in accordance with the purpose of the invention, as embodied and broadly described, a method according to the invention for registering a template image and a target image comprises several steps, including defining manifold landmark points in the template image and identifying points in the target image corresponding to the defined manifold landmark points. Once these points have been identified, the method includes the steps of computing a transform relating the defined manifold landmark points in the template image to corresponding points in the target image; fusing the first transform with a distance measure to determine a second transform relating all points within a region of interest in the target image to the corresponding points in the template image; and registering the template image with the target image using this second transform.
Therefore it is another embodiment of the present invention to provide a new framework for fusing the two separate image registration approaches in the large deformation setting. It is a further embodiment of the large deformation method to allow for anatomical and clinical experts to study geometric properties of imagery alone. This will allow for the use of landmark information for registering imagery. Yet another embodiment consistent with the present invention allows for registration based on image matching alone. Still another embodiment consistent with the present invention combines these transformations which are diffeomorphisms and can therefore be composed.
An embodiment consistent with the present invention provides an efficient descent solution through the Lagrangian temporal path space-time solution of the landmark matching problem for situations when there are small numbers of landmarks N less than  less than [xcexa9], reducing the xcexa9xc3x97T dimensional optimization to N, T-dimensional optimizations.
Another embodiment consistent with the present invention provides an efficient FFT based convolutional implementation of the image matching problem in Eulerian coordinates converting the optimization from real-valued functions on xcexa9xc3x97T to [T] optimizations on xcexa9.
To reformulate the classical small deformation solutions of the image registrations problems as large deformation solutions allowing for a precise algorithmic invention, proceeding from low-dimensional landmark information to high-dimensional volume information providing maps, which are one-to-one and onto, from which geometry of the image substructures may be studied. As the maps are diffeomorphisms, Riemannian lengths on curved surfaces, surface area, connected volume measures can all be computed guaranteed to be well defined because of the diffeomorphic nature of the maps.
It is also an embodiment consistent with the invention to introduce periodic and stationary properties into the covariance properties of the differential operator smoothing so that the solution involving inner-products can be implemented via convolutions and Fast Fourier transforms thus making the fusion solution computationally feasible in real time on serial computers
Both the foregoing general description and the following detailed description are exemplary and explanatory and are intended to provide further explanation of the invention as claimed.